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Abstract The Hard X-ray Modulation Telescope (HXMT) will perform an all-sky 
survey in hard X-ray band as well as deep imaging of a series of small sky regions. 
We expect various compact objects to be detected in these imaging observations. 
Point source detection performance of HXMT imaging observation depends not 
only on the instrument but also on its data analysis since images are reconstructed 
from HXMT observed data with numeric methods. Denoising technique plays an 
import part in HXMT imaging data analysis pipeline alongside with demodula¬ 
tion and source detection. In this paper we have implemented several methods 
for denoising HXMT data and evaluated the point source detection performances 
in terms of sensitivities and location accuracies. The results show that direct de¬ 
modulation with 1-fold cross correlation should be the default reconstruction and 
regularization methods, although both sensitivity and location accuracy could be 
further imporved by selecting and tuning numerical methods in data analysis of 
HXMT imaging observations. 

Key words: methods: data analysis — methods: numerical — techniques: image 
processing — instrumentation: high angular resolution 


1 INTRODUCTION 


1.1 Background 

Hard X-ray Modulation Telescope (HXMT) is a planned Chinese space X-ray telescope. The 
telescope will perform an all-sky survey in hard X-ray band (1 - 250 keV), a series of deep 
imaging observations at small sky regions and pointed observations. 

We expect a large number of X-ray sources, e.g., AGNs, to be detected in its all-sky survey. 
We also expect through a s eries of deep im aging observations at the Galactic plane X-ray 
transients can be detected!Li 2007 Lu 2012). Therefore the point source detection performance 


is one of our concerns on HXMT data analysis. 

Methods and corresponding sensitivities of pointed observation have been discussed by |Jin| 


et al. (2010). In imaging observations such as all-sky survey and deep imaging at small sky 


regions, a variety of data analysis aspects and methods are involved. First, images are com¬ 
puted instead of recorded directly by the optical instrument. Mapping as well as image 
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reconstruction methods are useful. The direct scientific data products from imaging observa¬ 
tions of HXMT are still scientific events, more specifically, X-ray photon arrival events. The 
attitude control system (ACS) of the spacecraft reports the attitudes periodically. We take these 
reported attitudes as nodes to perform certain interpolations to determine the instantaneous 
attitude for each scientific event. In this way, a set of parameters are assigned to each event, 
including the coordinates on the celestial sphere where the telescope is pointing at. Hence we 
call this process events mapping, where scientific events are mapped from time domain to the 
celestial sphere. The product of this process is refered to as the observed image, which implies 
the dimensionality of the data. Second, the exposure to a specific source is limited more 
strictly thus the signal-to-noise ratio (SNR) is tightly restricted. Hence the importance of reg¬ 
ularization methods becomes prominent. Finally, a picture is worth a thousand words. 
Various information can be extracted from an image by numerical methods. 

In this paper we investigate the point source detection performance of the imaging and 
detecting system synthesised from the telescope as well as diverse combinations of data analysis 
methods, especially the regularization methods. 

1.2 Modulation in HXMT imaging observation 

The PSF of HXMT HE, which is a composite telescope consisting of 18 collimated detectors, 
describes the response of the telescope to a point source when the telescope is pointing at the 
source as well as different locations around it. In other word, the PSF is a density function 
of the distribution of responses occur in different observation states, which is denoted by the 
instantaneous attitude of the telescope as a spacecraft. So the PSF takes the attitude of the 
telescope as its input. 

We use the proper Euler’s angles to describe the attitude of the telescope, i.e., </> and 9 
the longitude and latitude of the pointing, as well as ip denoting the rotation angle of the 
telescope around its own pointing axis, namely, the position angle. The modulation equation 
that corresponds to the imaging observation over the sphere surface is 

d{<p, e,iP)= IIp{<P, 9, ip, </>', 9')f{<p', 9') cos 9'd<P'A9', (1) 

n 

where f{(p', 9') is the object function (i.e., the image) defined in a compact subset of the sphere 
surface H, p{(p,9,ip,(p',9') is the modulation kernel function that relates the value of object 
function f{4>, 9') defined on a neighbourhood of the point (0', 9') of the subset H to the instan¬ 
taneous response of the telescope, while its status is {(p,9,ip). 

The modulation kernel function is determined by the point spread function (PSF) of the 
telescope. Suppose a unit point source is located at the zenith of the celestial sphere, i.e., the 
point (0, 0,1) in the corresponding Cartesian coordinate system. Fix the position angle ip, while 
slew the telescope across the polar cap, and assign responses of the telescope to the unit source 
to pixels on the celestial sphere. Then we obtain a map P{(p,9) represents the PSF on the 
celestial sphere with fixed rotation angle ip. 

The map is then projected to a tangent plane of the celestial sphere z = 1 hy gnomonic 
mapping, i.e., 

( u = cot 9 cos (p , . 

i , ( 2 ) 

( u = cot 9 sin (p ’ 

where u and v are local Cartesian coordinates on the tangent plane. Now we have 

Ptan)^, v) = Ptan(cot 9 COS 0, cot 9 sin 0) = P(0, 9), (3) 

i.e., the PSF dehned on the tangent plane. 
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To provide for data analysis we have two discrete forms of Ptan, 

, Pta.niu,v)dudv 


= 


Ila, , dwdr; 


( 4 ) 


where aij is the neighbourhood of the pixel (*,j) of the 2-D pixel grid, and the normalized 
form 

JL- . Pti,n{u,v)dudv 
JJ Pti,n{u,v)dudv 

Given the discrete image F, the detection area A and the duration of exposure on each 
pixel r, the observed data is 


D^t 


A 


max Pi 
ij 


{F * H) + ri, 


( 6 ) 


where rj, is the constant background count rate, * denotes the convolution, and H is the 
normalized PSF on the tangent plane. Eq. approximates modulations in HXMT imaging 
observations. 

Distortion occurs when projecting a set of points from (or to) a sphere surface to (or 
from) a plane. For example, distance between any two points, area of any continuous subset, 
angle between any two crossing lines (or tangents of curves) are altered non-uniformly. On the 
other hand, the rotation angle tp is not always fixed during HXMT imaging observations. Both 
distortions in image reconstruction from HXMT observed data and rotations during imaging 
observations have been discussed by Huo & Zhou (2013). Here for the sake of simplicity, we 
ignore them in this article. 

2 NUMERICAL METHODS 

2.1 Single point source detection performance estimation 

We estimate the single point source performance in terms of sensitivity and position accuracy 
through the following procedures. 

1. Determine flux threshold for point source detection. 


(a) 

(b) 

(c) 

(d) 


Simulate a frame of observed data contains only background counts. 

Run denoise program on the simulated data to try to increase the signal-to-noise ratio. 

Demodulate the denoised data. 


Run SExtractor, a source detection program, by Bertin & Arnouts (1996) on the de¬ 


modulated image to detect point sources and extract their intensities, coordinates and 
other parameters. At this point a catalog of point sources is compiled from the simu¬ 
lated data. Point sources detected here, i.e., from images demodulated from background 
data are false detections. 

(e) Repeat previous steps (from la to Id) and a series of catalogs are compiled. Draw a 


(f) 


histogram of flux of false detections that could possibly been detected from background 
counts given a specific condition of both observation and detection. 

Choose a cut from the histogram as the flux threshold so that a certain percent of 
the false detections are rejected and the rejection percentage is precise enough. The 
rejection percentage, e.g., 95% or 99.7% etc., reflects the significance of detections above 
the corresponding threshold. 

2. Estimate detection efficiency and position accuracy of a point source of specific flux intensity. 
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(a) 

(b) 

(c) 


,yr 


Simulate observed data contains a single point source of given flux intensity fm located 
j) in the model image. 

A catalog is compiled. 

Examine each detection is the catalog. Provided the i—th source in the catalog is de¬ 
tected at (xi,yi) in the demodulated image, the distance between the extracted source 
and the true point source 


at (a 

Perform steps from la to Id 


< 5 ,= 



-Xif + (jj 


m 



(7) 


as well as the flux of the extracted source fi are investigated to determine whether 
the i—th source is a true source or not. We define the score of the current catalog in 
detecting the single point source as 


dk = 


3 i : < A) A {fi > Fthres) 

otherwise 


( 8 ) 


where k is the index of the current catalog, A and Fthres are position accuracy and 
flux thresholds, therefore the score dk reveals whether the k—th. catalog contains the 
true source or not, in other word, through the previous steps (simulated observation, 
denoising, demodulation, source extraction and thresholding) whether we have detected 
the true source effectively. If we have, the outcome of these steps is counted as an 
effective detection of the true source, otherwise it is ineffective. 

(d) Repeat the previous steps (from 2a to 2c) N times and calculate the percentage of 


effective detections amongst all detections, namely. 


N 




k'l 


(9) 


which is defined as the detection efficiency. 

Let {xk,yk) be the position of the brightest source in the k—th catalog, the position 
accuracy is calculated as 


P = 



( 10 ) 


3. Find an flux intensity Fq.s so that the corresponding detection effeciency rj = 50%. The 
intensity Fq ^ marks the point source sensitivity of the detecting system synthesised from 
both the telescope in specific status and the data analysis program chain. 


2.2 Imaging and mapping 

2.2.1 Demodulation 


Direct demodulation (DD) method (Li & Wu 1994) is used to est imate the true image from 
observed data. Residual map calculated with CLEAN algorithm ( [Hogbom' 1974) is used as 
lower limit constraint in DD method. The skewness of the residual map is calculated in each 
CLEAN iteration and its minimum absolute value serves as the main stopping criterion of 
iterations. 
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2.2.2 Cross-correlation 


In contrary to detecting and extracting point sources from demodulated images, it is also feasible 
to do so from cross-correlated maps as long as the point sources are isolated with each other, 
compared to the FWHM of the PSF, since the position of a peak of the expected value of 
such a map coincides with the position of a source regardless of the PSF. Cross-correlating the 
observed data and the PSF yields the correlated map 

C = DirH, (11) 


where star denotes cross correlation, D and H denote the observed data and the PSF respec¬ 
tively. The peak of C that coincides with a point source of flux intensity F is 


C = 


T ■ A - F ^ 2 

maxi/„-^ 


( 12 ) 


according to Eq. On the other hand, the background variance of the correlated map is 




= T-rb-Y. 


(13) 


since r • rj, follows Poisson distribution. Hence the significance of the peak in term of numbers 
of cr is 


F-A-^r- E,,, HF 
maxi?* 


f-a-Vt 

\fFb 



(14) 


where T is total duration of exposure on the 2-D pixel grid, a/( i^tan) is the square root of the 
arithme mean of Ptan over the 2-D pixel grid, which is determined by the PSF as well as the 
range of the pixel grid only, provided the pixel grid is fine enough (see Eq. . 

Cross-correlation significance of isolated point source defined Eq. [T^ can be evaluated di¬ 
rectly, given only the flux intensity of the source, the background count rate, the duration of 
exposure, the detection area and the PSF. Hence it is determined by the object (i.e., the point 
source), the telescope and the status of observation thus effects from data analysis programs 
are minimized. 

We use the cross-correlation significance as a reference. For example, in our simulations an 
isolated point source of 1 mCrab flux has 2.42(T significance. 


2.3 Denoising 

2.3.1 Linear methods 


Gaussian smoothing is often used in digital image processing to suppress the noise at the cost 
of reduction in resolution. Trade-off between noise suppression and resolution conservation is 
adjusted through the standard deviation a of the Gaussian distribution serving as the smoothing 
kernel function. The best resolution of HXMT HE observed data is about 1.1°, limited by the 
FWHM of its narrow-field collimator. We set a to 28' so that the FWHM of the Gaussian kernel 
is also 1.1°. 

iV—fold cross correlation transform {N > 1) can be used in DD method to regularize the 
ill-posed problem, more specifically speaking, to ensure the convergence as well as stability of 


the solution! Li 2003). Here we put this technique in the denoising category. We have tested 


1—fold and 2—fold cross correlated DD methods respectively in this article. 
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2.3.2 Non-linear methods 


Non-local means denoising (Buades et al. 2005) is an edge-preserving non-linear denoising 


method. To increase its performance we have implemented this method with fast fourier trans¬ 
forms (FFTs). The pixel-wisely evaluation of the general Euclidean distance between the i—th 
pixel and other pixels of an image is replaced by 


D^ = {J2 ■Wk + {I^*W)-2[I^ {N, ■ IF)] 


(15) 


where Ni is neighbourhood of the *—th pixel, Ni^k is the k— pixel in the neighbourhood, I is 
the image, W is weight coefficient of the distance function. We use 7x7 pixels Gaussian kernel 
with standard deviation cr = 2 (in pixels) as the weight coefficient W in our simulations. We 


reduce the complexity of NLMeans method by computing convolutions in Eq. 15 with FFTs 
instead of using searching windows. ( jBuades et al. 2008). 

Median filter is another non-linear edge-preserving denoising method. This method is effec¬ 
tive in removing salt-and-pepper noise in digital images. In HXMT observed data such noise is 
incurred typically by missing data or charged particles of cosmic rays. We fixed the size of the 
filter at 2° x 2° (about 50 x 50 pixels) in HXMT HE data denoising. 

The last non-linear denoising method included in this article is the adaptive wavelet thresh¬ 
olding with multiresolution support (Starck & Murtagh 2006). The multiresolution support of 


a noisy image is a subset that contains significant coefficients only. So wavelet coefficients that 
dominated by noise are discarded. In this article we implemented the non-iterative algorithm. 
The 5 X 5 i ?3 spline wavelet is used for multiresolution decomposition. 


3 SIMULATION AND RESULTS 
3.1 In-orbit background simulation 

HXMT HE in-orbit background count rate rj, ranges from 147.6 counts/s to 210.7 counts/s 
(Li et al. 2009). We use a constant count rate 180 counts/s to simulate the average in-orbit 


background of HXMT HE in this article. 


3.2 Source energy spectrum and telescope detection efficiency 


We use the formula proposed by Massaro et al. (2000) together with parameters fitted by 
Jourdain & Roques (2009) 


= 3 87 X 75-1-79-0.1341ogio(E/2o) 


(16) 


to model energy spectra of Crab-like sources from 20 keV to 250 keV, where E is in keV and 
the flux F{E) is in photons/s/cm^/keV. 

The detector efficiency of HXMT HE is derived from its simulated energy response, as 
shown in Fig.[Tj Detection efficiency of HXMT HE to a Crab-like source is 67%. Count rate of 
HXMT HE corresponds to the 1 Crab intensity is 1112 counts/s, given the detection area of 
HXMT HE is about 5 100 cmF 


3.3 PSF and modulation 


We use the diagram shown in Fig. to simulate the PSF Ftan(it,u) on the tangent plane. 
We use the concentric average of P{(j),6), namely. 


s{0) = 


/:^p(0,0)d<(. 


27r 


(17) 
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Fig. 1 HXMT HE detection efficiency 


Simulated HXMT HE PSF 

Bo.g 

|o.8 




- 0.0 5 

- 0 . 10 

0.05 0.10 


S{0')d9' (18) 

to characterize the radial fade-out of the PSF and the concentration of the PSF respectively, 
as plotted in Fig. From Fig. we see that the FWHM of the simulated PSF is about 1.7° 
while 99.7% responses occur in a diameter of 11°. 

Despite the fact that the direct observed data is scientific events instead of 2-dimensional 
images, we start our simulation from simulated observed data in the form of images defined on 
2-dimensional cartesian pixel grid. We use a 22° x 22° model image for simulations. Given the 
diameter of the PSF, the central 11° x 11° region is fully modulated, that is, all contributions to 


- 0.10 - 0.05 0.00 


Fig. 2 Simulated HXMT HE PSF 


and the cumulative sum of S{0), 


C{0) = f 
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9 , in degrees 


Fig. 3 Radial distribution of the simulated HXMT HE PSF 


observed data in this region are from the model image only. The surrounding 33° x 33° region 
is partially modulated, i.e., part of the contributions to observed data is this region are from 
the model image. 

The average exposure per unit solid angle is 382 s/deg^ in HXMT half-year all-sky survey. 
The partially modulated region is discretized by x TV pixel grid, so r « 382 x 33^/jv^. The 
detection area of each HXMT HE detector is approximately 300 cm^ so the total area of all the 
17 HXMT HE detectors is 5 100 cm^. 


3.4 Results 


We have implemented several denoising methods (see Section 2.3). Fig.|^shows denoising results 
by these methods. 

We have simulated 5 000 frames of observed data that contains the in-orbit background 
counts only for each methods to estimate the corresponding flux thresholds by the method 
specified in Procedure]!] of Section 2.1 
thresholds, see Table JlJfor results. 


From false detections we have obtained 2tT and 3 (t flux 


Denoising 

2(7 thres., in 

mCrab 

3(7 thres., in 

mCrab 

w/o denoise 

0.520 

± 

0.001 

0.963 

± 

0.003 

1-fold CCT 

0.667 

± 

0.003 

1.192 

± 

0.010 

2-fold CCT 

0.919 

± 

0.014 

1.656 

± 

0.058 

Gaussian, a = 28' 

0.669 

± 

0.003 

1.133 

± 

0.008 

NLMeans 

0.845 

± 

0.004 

1.340 

± 

0.011 

Median filter, 2° 

0.709 

± 

0.004 

1.185 

± 

0.015 

Median filter, 4° 

0.665 

± 

0.005 

1.114 

± 

0.012 

Wavelet thres. 

0.0518 ±C 

.0002 

0.216 

± 

0.002 


Table 1 2cr and 3(7 thresholds of point source detection 


We have simulated 1000 frames of observd data that contains a Crab-like point source of 
1 mCrab, 1.25 mCrab, 1.5 mCrab, 1.75 mCrab, 2 mCrab, 2.5 mCrab, 3 mCrab, 4 mCrab, 
5 mCrab and 10 mCrab respectively, i.e., 10 000 frames of observed data in total. 
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Fig. 4 Denoising methods. From top to bottom, left to right: 1. Model image; 2. 
Observed data of 10 mCrab point source. 3. Gaussian smoothed data, a = 28'. 4. 
1-fold cross correlated data. 5. NLMeans filter denoised data. 6. 2° x 2° median filter 
denoised data. 7. 4° x 4° median filter denoised data. 8. 3cr wavelet thresholding 
denoised data, with spline wavelet transform. 


With these simulated data we have estimated location accuracies as well as detection effi¬ 
ciencies by the method described in Procedure of the following methods: 

1. DD without denoising, 

2. DD with 1-fold cross correlation, 

3. DD with 2-fold cross correlation, 

4. DD with Gaussian smoothing, where a = 28', 

5. DD with NLMeans filtering, the size of the filter is 7 x 7 and cr = 2 (both parameters are 
in pixels), 

6. DD with 2° x 2° median filter, 

7. DD with 4° x 4° median filter, and 
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8. DD with adaptive wavelet thresholding. 


Implementation details of the above methods are in Sect. 2.3 


Table shows the location accuracies on simulated data of 1 mCrab, 
and 10 mCrab point sources. 


2 mCrab, 5 mCrab 


1 mCrab 2 mCrab 5 mCrab 10 mCrab 


w/o denoise 

95 

± 

6 

39 

± 

2 

9.6 

± 

0.2 

5.2 

± 

0.1 

f-fold CCT 

104 

± 

3 

53 

± 

1 

20.5 

± 

0.4 

10.6 

± 

0.2 

2-fold CCT 

149 

± 

5 

90 

± 

2 

36.9 

± 

0.8 

17.2 

± 

0.4 

Gaussian, a — 28' 

90 

± 

4 

38 

± 

1 

11.1 

± 

0.2 

6.0 

± 

0.1 

NLMeans 

107 

± 

4 

52 

± 

1 

21.0 

± 

0.3 

12.2 

± 

0.2 

Median filter, 2° 

91 

± 

4 

43 

± 

1 

14.3 

± 

0.3 

7.7 

± 

0.1 

Median filter, 4° 

100 

± 

3 

57 

± 

1 

23.3 

± 

0.4 

14.8 

± 

0.3 

Wavelet thres. 

108 

± 

5 

48 

± 

2 

11.6 

± 

0.2 

5.9 

± 

0.1 


Table 2 Location accuracies, in arc minutes. 


Tablej^shows the detection efficiencies on simulated data of 1 mCrab, 1.5 mCrab, 2 mCrab, 
2.5 mCrab and 3 mCrab point sources. 

Although the RL iteration itself we employed in DD method is total-counts-conservative, 
i.e., the sum of counts of each pixel is conserved after the iteration ( Richardson[ |1972[ ), the 
regularizations including the background constrains as well as various denoising techniques 
are not necessarily counts-conservative. As a result, the absolute flux threshold (Table for 
rejecting false detections doesn’t reflect the sensitivity directly but the detection efficiency 
(Table [|) does. 

Errors of flux thresholds, location accuracies and detection efficiencies in Table Table 
and Table are calculated by bootstrapping. For example, following Procedure a set of 5 000 
frames of demodulated images are obtained, from which false detections are calculated and a 
histogram is plotted, where both the 2a and 3a thresholds are determined. Now let’s generate 
a new set with the same volume by resampling from the original set with replacement in order 
to calculate both the thresholds again. We repeat this resampling process until we get enough 
thresholds calculated to estimate their standard deviations. In Table Table and Table 
each of the errors is a standard deviation that calculated from 1000 resampled sets. 


f mCrab 1.5 mCrab 2 mCrab 2.5 mCrab 3 mCrab 


w/o denoise 

29 

± 

2 

53 

± 

2 

77 

± 

1 

94 

± 

1 

98 

± 

0 

l-fold CCT 

41 

± 

2 

71 

± 

1 

92 

± 

1 

99 

± 

0 

100 

± 

0 

2-fold CCT 

27 

± 

2 

56 

± 

2 

79 

± 

2 

95 

± 

1 

99 

± 

0 

Gaussian, a — 28' 

35 

± 

2 

66 

± 

2 

87 

± 

1 

98 

± 

0 

100 

± 

0 

NLMeans 

36 

± 

2 

68 

± 

1 

91 

± 

1 

100 

± 

0 

100 

± 

0 

Median filter, 2° 

37 

± 

2 

64 

± 

2 

86 

± 

1 

97 

± 

1 

100 

± 

0 

Median filter, 4° 

39 

± 

2 

69 

± 

1 

89 

± 

1 

98 

± 

0 

100 

± 

0 

Wavelet thres. 

27 

± 

1 

51 

± 

2 

79 

± 

1 

95 

± 

1 

100 

± 

0 


Table 3 Detection efficiencies, in percentages. 


A comprehensive summary of all the tested methods on all simulated data is shown in Fig. 

El 

4 CONCLUSION 

According to the results from our tests no denoising method shows significant advantage over 
1-fold cross correlated DD method in single point source detection efficiency. Therefore it’s 
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Fig. 5 Single source location accuracy and detection efficiency. MF stands for median 
filter while WT for wavelet thresholding. 


suggested that 1-fold cross correlation should be the default regularization method for single 
point source detection in HXMT imaging data analysis. 

In the other hand, the location accuracy can be improved with alternative denoising meth¬ 
ods, such as median filter, wavelet thresholding, Gaussian smoothing with little kernel, or with¬ 
out denoising, according to the results in this work. 

This article is focused on the single point source detection of HXMT imaging data analy¬ 
sis, where other interesting topics can not be all covered. Although the alternative denoising 
methods have been out-performed more or less by the default 1-fold cross correlation in their 
contributions to the detection efhciency, they have shown certain advantages in location accu- 
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racy, and these features are promising for locating bright transients, multiple sources resolving, 
and so on. 
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